% Routine that searches over parameter values by annealing algorithm.

function [mincriterion Optparams optmom results NSE] = annealing(init, span, nK, nL, nA, nF, extraA, Lmode, FC, momlist, Spec)

%          [(1)    (2)    (3)    (4)    (5)       (6)     (7)    (8)    (9)    (10)   (11)  (12)  (13)]     
% params = [ bK     bL   alphaK alphaL alphaKL resaleloss rho   sigma   sigdK  sigdL   ro   meK   meL ]
upcap    = [7.50   7.50   0.25   0.25   0.25      1.00    0.99   0.80   0.0006 0.0006 1.00  1.00  1.00];
lowcap   = [0.00   0.00   0.00   0.00  -0.25      0.00    0.00   0.01   0.0000 0.0000 1.00  0.00  0.00];
% Initial guess:
  params  = init;
  % Resulting criterion evaluation:
  [criterion simmom] = simulation(params, nK, nL, nA, nF, extraA, Lmode, FC, 'Search', momlist);
  if isnan(criterion) == 1;
     criterion1 = 1e7;
  else
     criterion1 = criterion;
  end
  simmom1 = simmom;
  params1 = params;
  guess = 0.0; 
  scalar = 0.50;

  % Start the annealing loop:
  for ii = 1:1000;
    jj = ii^(1/(1+sqrt(sum(span>0))));
    RandStream.setDefaultStream(RandStream('mt19937ar','seed',sum(100*clock)));
    z    = span.*(1/jj).*randn(1,length(params));
    params=min(max(params+z+guess,lowcap),upcap);
    if params(5) < 0;
        params(5) = max(params(5),-min(params(3),params(4))); % Make shure that alphaKL does not take a negative value that exceeds the size of alphaK or alphaL. This can be implemented and made more flexible in simulation.m later on. 
    end
    
    [criterion, simmom, ~, ~, ~, ~, ~, ~, ~, ~, ~, ~, ~, ~, ~] = simulation(params, nK, nL, nA, nF, extraA, Lmode, FC, 'Search', momlist);
    
    % Show the results after the loop:
    
    fprintf('---------------------------------------- \n')
    fprintf('---------------------------------------- \n')
    fprintf('Results: Iteration # %3.0f  \n',ii)
    fprintf('---------------------------------------- \n')
    fprintf('Parameter   Current     Best \n')
    fprintf('---------------------------------------- \n')
    fprintf('bK         = %2.4f    %2.4f   \n', params(1), params1(1))
    fprintf('bL         = %2.4f    %2.4f   \n', params(2), params1(2))
    fprintf('alphaK     = %2.4f    %2.4f   \n', params(3), params1(3))
    fprintf('alphaL     = %2.4f    %2.4f   \n', params(4), params1(4))
    fprintf('alphaKL    = %2.4f    %2.4f   \n', params(5), params1(5))
    fprintf('RSloss     = %2.4f    %2.4f   \n', params(6), params1(6))
    fprintf('rho        = %2.4f    %2.4f   \n', params(7), params1(7))
    fprintf('sigma      = %2.4f    %2.4f   \n', params(8), params1(8))
    fprintf('sigdK      = %2.4f    %2.4f   \n', params(9), params1(9))
    fprintf('sigdL      = %2.4f    %2.4f   \n', params(10), params1(10))
    fprintf('ro         = %2.4f    %2.4f   \n', params(11), params1(11))
    fprintf('meK        = %2.4f    %2.4f   \n', params(12), params1(12))
    fprintf('meL        = %2.4f    %2.4f   \n', params(13), params1(13))
    fprintf('---------------------------------------- \n')
    fprintf('nK         = %3.0f  \n', nK)
    fprintf('nL         = %3.0f  \n', nL)
    fprintf('nA         = %3.0f  \n', nA)
    fprintf('nF         = %3.0f  \n', nF)
    fprintf('extraA     = %3.0f  \n', extraA)
    fprintf('---------------------------------------- \n')
    fprintf('---------------------------------------- \n')
    fprintf('Criterion             = %6.3f \n', criterion)
    fprintf('Best criterion        = %6.3f \n', criterion1)
    fprintf('Inv(Chi2(0.95)) (df)  = %6.3f (%1.0f) \n', icdf('chi2',0.95,length(momlist)-sum(span>0)),length(momlist) - sum(span>0))
    fprintf('---------------------------------------- \n')
    fprintf('---------------------------------------- \n \n \n')
    
    if isnan(criterion)==1;
        params = params1;
        guess=scalar*(params1-params);
    else
        if criterion < criterion1;        % If evaluation after current guess is better than best guess so far.
            guess=scalar*(params-params1);  % Determines the direction of next guess.
            params1=params;                 % Updates the preferred set of parameters so far.
            criterion1=criterion;           % Updates the best evaluation with result after best guess so far. 
            simmom1=simmom;                 % Updates the best fitted simulated moments so far.
        else                              % If evaluation after current guess is not better than best guess so far.
            guess=scalar*(params1-params);  % Determines the direction of next guess.
            params=params1;                 % Resets the choice of paramters to best guess so far. 
        end
    end
  end    
    
  % After the annealing loop has stopped searching, we can go ahead and calculate numerical standard errors. We call sterrors.m to do this.
  Optparams = params1;
  save('workspace.mat','Optparams')
  % NSE = zeros(12,1);
  [results NSE] = sterrors(params1, criterion1, simmom1, nK, nL, nA, nF, extraA, Lmode, FC, momlist, Spec);
  optmom = simmom1;
  mincriterion = criterion1;
